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A  Project  Abstract 


To  meet  the  military’s  increased  need  for  rapidly  deployable  communication  solutions,  wireless 
networks  are  becoming  increasingly  common  within  tactical  environments.  Frequency  hopping 
(FH)  is  widely  used  for  radio  transmission  in  such  networks  due  to  its  low  probability  of  de¬ 
tection/interception.  With  the  increasing  deployment,  multiple  networks  occupying  overlapping 
frequency  bands  are  likely  to  coexist  in  a  physical  environment,  especially  in  tactical  operations, 
emergency  situations,  or  dense-populated  areas.  Consequently  co-channel  interference  due  to  fre¬ 
quency  collisions  can  become  a  major  performance  limiting  factor.  In  this  project,  we  developed  a 
novel  interference  mitigation  technique  based  on  multidimensional  frequency  estimation  coupled 
with  the  expectation-maximization  principle,  which  effectively  resolves  collisions  among  non- 
collaborative  networks,  thus  enabling  robust  coexistence  of  multiple  FH  networks.  To  deal  with 
possible  receiver-transmitter  mismatch,  we  also  designed  a  low-complexity  model  order  variation 
detection  method.  Novel  multidimensional  frequency  estimation  algorithms  are  also  investigated. 
The  significance  of  this  project  in  basic  research  lies  in  innovative  schemes  for  collision  resolution 
enabling  interference-robust  FH  networking. 


B  Technical  Report 

B.l  Problem  Statement 

To  meet  the  military’s  growing  need  for  rapidly  deployable  communication  solutions,  wireless  net¬ 
works  are  becoming  increasingly  common  in  tactical  environments.  For  example,  military  units 
(e.g.,  soldiers  and  tanks),  equipped  with  wireless  devices,  could  form  multiple  networks  in  tactical 
operations  when  they  roam  in  a  battlefield.  Frequency  hopping  spread  spectrum  (FHSS)  is  widely 
used  for  radio  transmission  in  such  networks,  due  to  its  low  probability  of  detection/interception 
[1].  For  example,  FHSS  is  adopted  in  the  single  channel  ground-airborne  radio  system  (SINC- 
GARS),  which  is  the  current  standard  army  combat  net  radio.  Recently  FHSS  has  also  been 
adopted  in  commercial  applications  such  as  Bluetooth  [2]. 

Because  of  their  increasing  deployment,  multiple  (homogeneous  and/or  heterogeneous)  wire¬ 
less  networks  with  overlapping  frequency  bands  are  likely  to  coexist  in  a  physical  environment, 
especially  in  tactical  operations,  emergency  situations,  or  dense-populated  areas.  Consequently 
co-channel  interference  due  to  frequency  collisions  can  become  a  major  performance  limiting  fac¬ 
tor  [3-5].  When  collisions  occur,  network  throughput  decreases  and  delay  can  become  excessive 
due  to  retransmissions.  Theoretical  analysis  has  shown  that  the  packet  error  rate  (PER)  of  one 
Bluetooth  piconet  due  to  collisions  can  be  up  to  10%  if  seven  piconets  coexist  [6],  and  a  Bluetooth 
receiver  may  experience  up  to  27%  packet  loss  for  data  traffic  in  the  presence  of  interference  from 
an  IEEE  802.11b  based  wireless  local  area  network  (WLAN)  [7]. 

Recently  the  coexistence  issue  has  gained  increasing  attention  [8-14].  However,  to  date  most 
coexistence  schemes  are  designed  for  simultaneous  functionality  of  a  Bluetooth  piconet  and  an 
802.11b  WLAN.  The  latter  is  a  direct  sequence  spread  spectrum  (DSSS)  network  that  occupies  a 
fixed  frequency  band  of  22  MHz.  In  this  project,  we  study  the  problem  of  coexistence  of  multi- 


1 


pie  FH  networks,  where  the  interfering  frequency  channels  are  constantly  changing  and  the  hop 
sequence  of  one  network  is  unknown  to  another.  We  develop  a  physical  layer  method  to  mitigate 
co-channel  interference  and  enable  robust  coexistence  of  multiple  FH  networks. 

B.2  Prior  Work 

Most  coexistence  methods  are  developed  for  simultaneous  function  of  a  Bluetooth  piconet  and 
a  WLAN,  based  on  various  MAC  layer  mechanisms  that  enable  the  Bluetooth  devices  to  avoid 
hopping  onto  the  frequency  band  occupied  by  the  WLAN.  Other  methods  mitigate  co-channel 
interference  with  physical  layer  approaches. 

B.2.1  Collision  Avoidance  for  Coexistence  of  Bluetooth  and  WLAN 

The  collision  avoidance  techniques  can  be  classified  as  collaborative  and  non-collaborative  ones. 
For  collaborative  cases,  attractive  data  transmission  rates  and  throughput  can  be  achieved  by  using 
a  communication  link  between  the  Bluetooth  and  WLAN  when  they  are  embedded  on  the  same 
device  [15]  or  by  coordinating  the  hop  frequencies  of  the  co-located  Bluetooth  devices  [8]. 

Non-collaborative  methods  do  not  require  direct  communication  between  the  two  networks, 
and  they  usually  rely  on  monitoring  the  channel  to  detect  interference  and  estimate  traffic.  For  ex¬ 
ample,  power  control  is  employed  based  on  PER  [9]  or  the  received  signal  strength  [10]  to  sustain 
the  quality  of  service  for  a  Bluetooth  link.  To  avoid  hopping  onto  preoccupied  frequency  chan¬ 
nels,  adaptive  frequency  hopping  (AFH)  modifies  the  Bluetooth  frequency  hopping  sequence,  and 
Bluetooth  interference  aware  scheduling  (BIAS)  strategy  postpones  the  transmission  [11],  both  de¬ 
tecting  preoccupied  frequency  bands  by  monitoring  PERs  on  all  channels.  The  overlap  avoidance 
(OLA)  schemes  proposed  in  [12]  are  based  on  packet  scheduling  and  traffic  control.  There  are 
also  hybrid  methods  that  combine  AFH  and  Bluetooth  carrier  sense  (BCS)  [13]  or  combine  power 
control,  listen-before-talk  (LBT)  and  AFH  [14]  to  achieve  improved  performance. 

The  center  control  mechanism  needed  for  collaborative  schemes  [8, 15]  confines  their  appli¬ 
cations  to  certain  situations.  Power  control  methods  [9, 10]  depend  on  the  accuracy  of  channel 
sensing  and  can  not  provide  much  improvement  if  the  Bluetooth  device  is  very  close  to  the  inter¬ 
fering  device.  Carrier  sensing  based  schemes  inevitably  suffers  from  the  hidden  terminal  prob¬ 
lem  [13, 16].  Approaches  based  on  scheduling  such  as  BIAS  [11],  OLA  [12]  and  master  delay 
MAC  scheduling  (MDMS)  [17]  cause  delay  in  the  transmission,  hence  they  may  not  be  bandwidth 
efficient.  AFH  [11]  and  interference  source  oriented  AFH  (ISOAFH)  [17]  are  effective  in  dealing 
with  WLAN  interference,  but  not  applicable  for  multiple  co-located  FH  networks.  The  perfor¬ 
mance  of  AFH  is  also  dependent  on  the  update  rate  of  the  frequency  classification  to  track  the 
channel  dynamics  [11].  A  hybrid  method  of  power  control,  LBT  and  AFH  proposed  in  [14]  can 
achieve  better  performance  but  will  add  complexity  to  the  application. 

In  summary,  most  non-collaborative  interference  detection  and  collision  avoidance  schemes  are 
developed  for  coexistence  of  Bluetooth  and  WLAN,  and  they  are  not  applicable  to  the  coexistence 
of  multiple  FH  networks,  because  the  frequency  channels  are  constantly  changing  and  the  hop 
sequence  of  one  network  is  not  known  to  another. 
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B.2.2  Co-channel  Interference  Mitigation  in  FH  Systems 

Instead  of  avoiding  the  collision  by  scheduling,  many  efforts  have  been  made  on  interference  sup¬ 
pression  and  mitigation  on  physical  layer.  The  problem  is  challenging  because  when  multiple  FH 
wireless  networks  coexist  in  noncooperative  scenarios,  not  only  the  hop  sequences,  hop  timing/rate 
and  other  parameters  such  as  hop  bandwidth  and  bin-width  are  unknown  to  each  other,  but  also 
there  maybe  model  order  variations  even  if  the  number  of  active  emitters  remains  the  same.  If  the 
receiver’s  bandwidth  is  not  matched  to  the  hop  bandwidth  of  the  emitters,  the  FH  signals  may  hop 
in  and  out  of  the  observation  frequency  band  of  the  receiver. 

Considerable  work  has  been  done  on  blind  interference  suppression  for  FHSS  systems  using 
an  antenna  array  [18,  19].  These  approaches  aim  at  interference  suppression  rather  than  joint 
multiuser  detection  and  hop  timing  estimation,  and  most  of  them  require  at  least  the  knowledge  of 
the  hop  sequence  of  the  user  of  interest,  and  their  interference  nulling  capability  is  bounded  by  the 
degrees  of  freedom  in  the  adaptive  array.  Multiuser  detection  for  FH  systems  have  been  considered 
in  [20,21].  They  assume  that  the  hop  sequences  of  are  known  to  the  receiver,  and  thus  they  are 
not  applicable  to  noncooperative  FH  emitters.  In  [22]  the  estimation  of  the  hop  sequence  for  a 
single  user  is  discussed,  while  the  remaining  users  are  treated  as  white  Gaussian  interference.  The 
approach  of  [22]  is  conceptually  simple,  but  assumes  perfect  channel  knowledge. 

Without  assuming  the  knowledge  of  hop  sequences,  several  (semi-)blind  methods  have  been 
proposed  for  hop  timing  and  frequency  estimation.  For  example,  with  known  hop  rate  and  fre¬ 
quency  bins,  channelized  receivers  have  been  proposed  for  semi-blind  hop  timing  estimation  [23, 
24].  However,  the  performance  of  those  receivers  degrades  rapidly  if  the  channelization  is  im¬ 
perfect,  for  example,  when  there  is  unknown  bandwidth  mismatch  or  carrier  frequency  offset. 
In  [25,  26],  hop  timing  and  frequency  estimators  based  on  the  principle  of  dynamic  program¬ 
ming  (DP)  were  developed  for  blind  tracking  of  multiple  FH  signals,  using  either  a  decoupled  ap¬ 
proach  [25]  or  a  joint  approach  [26].  These  methods  neither  assume  knowledge  of  hop  sequences 
or  hop  timing,  nor  rely  on  channelization,  and  hence  are  robust  to  frequency  offset.  However, 
they  require  accurate  model  order  information,  and  thus  can  not  handle  bandwidth  mismatch.  In 
addition,  the  complexity  of  DP  is  too  high  to  be  feasible  for  practical  implementation.  An  iterative 
maximum  likelihood  (ML)  algorithm  is  developed  in  [27]  to  estimate  the  hop  timing  and  frequen¬ 
cies  with  low  complexity,  but  it  can  only  deal  with  the  single-user  two-hop  case  (one  hop  timing 
instant  to  be  estimated). 

B.3  Summary  of  Major  Results 

In  this  project,  we  develop  a  physical  layer  method  for  co-channel  interference  mitigation  and 
robust  coexistence  of  multiple  FH  networks.  The  main  results  are  summarized  as  follows. 

1.  We  develop  an  expectation-maximization  (EM)  algorithm  combined  with  a  2-D  frequency 
estimator  for  hop  timing  and  frequency  estimation  of  multiple  FH  signals  with  transmitter- 
receiver  bandwidth  mismatch.  The  algorithm  resolves  frequency  collision  without  retrans¬ 
mission.  A  software  testbed  is  developed  by  applying  the  algorithm  in  Bluetooth  piconets. 
Simulation  results  demonstrate  the  effectiveness  of  the  proposed  method. 
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2.  When  there  is  transmitter-receiver  bandwidth  mismatch,  FH  signals  may  hop  in  and  out  of 
the  observation  band  of  the  receiver,  which  results  in  model  order  variation.  We  design  a  low 
complexity  approach  for  model  order  variation  detection  based  on  the  trace  of  the  covariance 
matrix  of  the  received  signals.  This  approach  outperforms  the  sliding-window  singular  value 
decomposition  (SVD)  or  minimum  description  length  (MDL)  methods. 

3.  Multidimensional  frequency  estimation  plays  an  important  role  in  collision  resolution  of 
multiple  FH  signals.  We  optimize  eigenvector-based  multidimensional  frequency  estimation 
by  minimizing  the  estimation  error  variance. 

In  the  following,  upper  (lower)  bold  face  letters  are  used  for  matrices  (column  vectors).  A*, 
Ar,  Ah ,  and  A'  denote  the  conjugate,  transpose,  Hermitian  transpose,  and  pseudo-inverse  of 

A,  respectively.  We  use  (8)  for  the  Kronecker  product,  and  0  for  the  Khatri-Rao  (column-wise 
Kronecker)  product.  We  also  use  Ip  for  a  p  x  p  identity  matrix,  0 mxn  for  an  M  x  N  zero  matrix, 
D  (a)  for  a  diagonal  matrix  with  a  as  its  diagonal,  A(m)  for  a  sub-matrix  of  A  formed  by  its  first 
m  rows,  and  ||  ■  ||  for  the  l2  norm. 

B. 4  Collision  Resolution  Using  an  Array  with  Bandwidth  Mismatch 

Relying  on  the  principle  of  expectation-maximization  and  efficient  2-D  frequency  estimation,  we 
develop  an  EM  algorithm  to  jointly  estimate  the  hop  timing  and  frequencies  of  multiple  frequency 
signals  without  the  knowledge  of  their  hop  patterns,  and  it  remains  operational  in  the  presence  of 
identical  frequencies  and  bandwidth  mismatch.  The  idea  behind  the  collision  resolution  algorithm 
is  that  if  a  hop-free  dataset  were  available,  one  could  model  the  collided  data  packets  as  a  mixture  of 
(modulated)  complex  exponentials.  Cast  in  proper  matrix  form,  such  a  signal  has  a  Vandermonde 
structure  in  the  time  domain.  In  addition,  the  use  of  a  uniformly  spaced  linear  array,  further  induces 
Vandermonde  structure  in  the  spatial  domain.  We  exploit  the  2-D  Vandermonde  structures  and  use 
a  2-D  frequency  estimation  algorithm  that  draws  upon  the  rich  identifiability  and  near-optimality 
results  in  [28]  to  recover  the  hop  frequencies.  Hop  timing  estimates  are  obtained  by  coupling  the 
2-D  frequency  estimation  with  the  expectation  maximization  algorithm,  thus  collided  packets  can 
be  recovered  without  retransmission. 

B.4.1  Data  Model 

Suppose  at  time  t  G  [tk,tk+ i),  an  M -element  uniform  linear  array  (ULA)  receives  a  total  of  ci/, 
effective  signals.  Each  far  field  FH  signal  is  from  a  nominal  direction-of-arrival  (DOA)  with  neg¬ 
ligible  angle  spread.  If  there  exists  transmitter-receiver  bandwidth  mismatch,  the  FH  signals  hop 
in  and  out  of  the  observation  band  of  the  receiver,  which  results  in  model  order  variations.  Here 
model  order  variation  refers  to  the  change  of  d fc.  The  M  x  1  baseband  received  signal  vector  at  the 
array  output  can  be  written  as 


=  E  4“  ^(^)?  ^  ^k+ 1?  (1) 

2=1 


4 


Frequency 


Figure  1 :  An  illustration  of  two  FH  signals 

where  tk  is  the  A  -th  system- wide  hop  instant  (hop  timing),  and  d /,  is  the  number  of  effective  signals 
during  the  time  period  delimited  in  tk  <  t  <  tk+1.  Here  a(9itk)  is  the  array  steering  vector,  which 
can  be  represented  as  a(6itk)  =  [l  •  •  ■  0fIk~1]T,  and  6hk  =  eJ27rAsm(“»,fc),  where  aijk  is  the  DOA 

of  the  i-th  narrow-band  signal  during  the  A-th  time  period,  and  A  is  the  array  baseline  separation 
in  terms  of  wavelength.  In  (1),  6^k  denotes  the  complex  path  loss  for  the  ?-th  signal  that  collects 
the  overall  attenuation  and  propagation  phase  shift.  v(t)  is  the  complex  white  Gaussian  noise 
vector.  For  simplicity  of  exposition,  the  transmitted  FH  signal  is  assumed  to  be  FSK  modulated, 
though  the  algorithm  can  also  accommodate  other  modulations  such  as  phase  shift  keying  (PSK)  or 
quadrature  amplitude  modulation  (QAM).  The  i-th  transmitted  FH  signal  is  s,j,(t)  =  The 

power  of  Sijk(t)  is  absorbed  into  3r±.  Here  the  carrier  shifts  due  to  hopping  or  symbol  modulation 
are  treated  as  conceptually  equivalent,  albeit  of  different  magnitudes. 

An  illustration  of  two  FH  signals  is  shown  in  Fig.  1,  where  tk  denotes  a  system-wide  hop 
instant.  The  two  FH  signals  are  asynchronous  and  have  different  hop  rates.  The  received  signal 
in  (1)  is  hop-free  between  any  two  adjacent  system-wide  hop  instants.  Signal  1  hops  out  of  the 
observation  band  at  tk  and  hops  back  in  at  tk+ 1,  so  the  number  of  the  effective  signals  in  (1)  during 
tk  <t  <  tk. |_i  is  dk  —  1,  while  during  other  hop-free  segments  the  model  order  is  2. 

After  sampling  the  signals  in  (1)  with  a  normalized  sampling  period,  the  discrete-time  equiva¬ 
lent  model  becomes 


dk 

xn  =  a(Oi,k)Pi,kSi,k(n)  +  vn ,  nk  <  n<  nk+i,  (2) 

i=  1 

where  n  —  0,  •  •  •  ,  N  —  1,  and  N  is  the  total  number  of  snapshots;  nk  is  the  A-th  hop  instant,  where 
A  =  0, . . . ,  K ,  and  K  is  the  total  number  of  hops.  By  default  we  let  no  =  0  and  rix  =  N.  The 
n-th  snapshot  of  the  z-th  transmitted  signal  is  sitk(n )  =  eju}i’k'n,  nk  <n<  nk+i,  where  ay/,  is  the 
frequency  of  the  z-th  signal  during  the  time  period  from  nk  to  nk+ ,  —  1.  Since  nk  is  a  system-wide 
hop  instant,  it  is  possible  that  u>ijk  =  ut^-i  for  some  i.  If  a’,./,  =  ay./,,  a  collision  occurs.  The 
objective  here  is  to  estimate  {nk}  and  hop  frequencies  for  a  given  set  of  observations  {xn}^=(' , 
with  unknown  DOAs,  hop  sequences,  amplitudes,  and  number  of  effective  signals.  After  these 
parameters  are  obtained,  collisions  are  then  resolved. 
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B.4.2  Known  Hop  Timing 


If  hop  timing  is  known,  i.e.,  nk  s  are  known,  then  the  received  data  set  can  be  split  into  K  hop-free 
data  subsets.  Between  any  two  consecutive  hop  instants,  e.g.,  nk  and  nk+ 1,  there  are  dk  temporal 
frequencies  involved.  During  such  a  time  segment,  the  k- th  received  data  subset  can  be  written  as 

Xk  =  AkBkS k  +  V  fc,  (3) 

where  Xk  =  [  xnk  xnk+\  •  •  ■  xnk+1_i  ] ,  and  Ak,  Bk,  Sk,  and  Vk  are  the  corresponding  matrices 
for  antenna  response,  channel  attenuation,  transmitted  signals,  and  noise.  Note  that  the  number 
of  effective  signal  is  dk  (i.e.,  model  order),  which  can  be  estimated  from  X k  by  an  appropriate 
source  enumeration  method  because  it  is  fixed  in  this  subset.  A  source  enumeration  method  can  be 
based  on  rank  detection  criteria  (e.g.  SVD  [29])  or  information  theoretic  criteria  (e.g.,  the  Akaike 
information  criterion  (AIC)  [30],  MDL  [31,32]).  SVD  requires  the  knowledge  of  noise  variance  for 
threshold  setting  in  model  order  selection  [29].  AIC  tends  to  overestimate  the  number  of  signals  in 
the  large  sample  limit,  while  MDL  is  shown  to  be  a  consistent  estimator  [31].  Therefore  we  choose 
MDL  to  estimate  dk  from  Xk,  which  is  to  calculate 

(  i  \  (M-m)(nk+1-nk) 

V\M  \M~m  \  1 

' — C  +  2m(2M  ~  m)  log(™fc+i  -  nk ),  (4) 

M—m  X/i=rn+\  *  J 

where  Ai  >  A2  •  •  •  >  A M  are  the  eigenvalues  of  the  sample  covariance  matrix  of  X k.  The  number 
of  signals  dk  is  determined  as  the  value  of  m  G  (0, 1, . . . ,  M  —  1}  for  which  MDL(m)  in  (4)  is 
minimized.  If  dk  >  M,  this  method  does  not  work.  However,  exploiting  the  inherent  structure  of 
the  data  matrix  Xk,we  can  use  data  smoothing  to  expand  the  size  the  X k  while  keeping  its  model 
order  unchanged.  Then,  the  MDL  algorithm  can  be  applied  to  the  expanded  data  matrix  to  estimate 
the  model  order  even  if  dk  >  M  (see  [33]  for  details). 

After  dk  is  obtained,  the  estimation  of  DOAs  and  frequencies  from  X  k  in  (3)  can  be  viewed  as 
a  2-D  frequency  estimation  problem  because  both  Ak  and  Sk  are  Vandermonde  matrices.  There 
are  dk  frequency  components  along  each  of  the  spatial  and  temporal  dimensions  of  the  2-D  fre¬ 
quency  mixture  Xk.  Various  subspace  methods  can  be  used  to  estimate  di>k,  A, a-,  and  uj^k,  for 
i  =  1, . . .  ,dk,  from  Xk.  For  example,  eigenvalue-based  algorithms  such  as  estimation  of  signal 
parameter  via  rotational  invariance  technique  (ESPRIT)  [34],  unitary  ESPRIT  [35],  joint  angle- 
frequency  estimation  (JAFE)  [36],  matrix  enhancement  and  matrix  pencil  (MEMP)  [37],  and 
eigenvector-based  algorithm  such  as  multidimensional  folding  (MDF)  [28],  are  all  applicable.  We 
choose  the  MDF  algorithm  here  because  MDF  achieves  the  most  relaxed  identifiability  bound  and 
is  shown  to  outperform  ESPRIT-like  algorithms  such  as  JAFE  and  MEMP,  and  do  not  require  an 
extra  frequency  pairing  step  (see  [28,38]  for  more  detail).  The  description  of  the  MDF  algorithm 
can  be  found  in  [28, 38]  and  is  omitted  here. 

B.4.3  Hop  Timing  Estimation 

Both  the  model  order  estimation  by  the  MDL  algorithm  and  2-D  frequency  estimation  by  the  MDF 
algorithm  need  to  work  with  a  hop-free  data  subset.  However,  here  the  hop  timing  is  unknown.  In 
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the  following,  we  develop  an  EM  algorithm  that  yields  joint  estimates  of  hop  timing  and  frequen¬ 
cies.  It  is  shown  in  [39]  that  the  EM  algorithm  can  be  written  as 


E-step:  for  1  <  k  <  K  —  1,  compute 


(p) 

nV  =  arg  mm 

nk 


-  (sf) 


+ 


V  A  (p)  r(p)  (  C (p) 

1  “  A-k-l&k-l  1 


where  nk  G 

M-step:  compute 


Ap) 


</>(,J+1)  =  arg  min  ( 

f  K-l 

^  .  A,  -  AkBkS 

l  fc=0 

||2f 

\\f 

f 

(5) 


(6) 


As  discussed  in  Sec.  B.4.2,  when  n(p>  is  given,  the  model  order  can  be  obtained  using  the  MDL 
algorithm  and  other  parameters  of  c/)l'p+1'1  can  be  obtained  by  applying  the  MDF  algorithm  to  the 
corresponding  X /.\s  determined  by  n(pK  It  can  be  observed  from  (5)  and  (6)  that  the  proposed  EM 
approach  is  actually  a  sort  of  decoupled  ML  algorithm,  since  n  and  <fi  are  estimated  in  the  E-step 
and  M-step  respectively  in  each  iteration. 

In  summary,  given  a  received  data  block  with  model  order  variation  and  multiple  hops,  the 
EM  algorithm  first  takes  a  randomly  generated  or  pre-determined  n([y  as  the  initialization  for 
the  EM  algorithm.  Then,  it  iterates  the  following  two  steps  until  convergence:  the  M-step,  Eqn. 
(6),  provides  a  new  estimate  of  model  orders  and  hop  frequencies,  which  are  accomplished  by 
applying  first  the  MDL  algorithm  and  then  the  MDF  algorithm  to  data  segments  according  to 
the  updated  assignment;  the  E-step,  Eqn.  (5),  involves  assigning  signal  segments  to  the  current 
estimated  parameters  that  fit  them  best.  Upon  convergence,  frequencies  and  complex  amplitudes 
of  different  segments  pertaining  to  a  particular  signal  can  be  associated  via  their  corresponding 
DOA  parameters,  since  for  a  single  data  subset,  frequency  and  DOA  parameters  pertaining  to 
one  signal  are  paired  up  automatically  by  the  MDF  algorithm.  Therefore,  joint  hop  timing  and 
frequency  estimation  for  multiple  FH  signals  is  achieved  and  collisions  are  resolved.  Furthermore, 
to  improve  the  performance,  in  [39]  we  also  develop  a  low  complexity  initialization  based  on  the 
spectrogram  of  the  received  data  for  the  EM  algorithm. 


B.4.4  Simulation  Results 

In  this  section,  we  first  compare  the  EM  algorithm  with  the  iterative  ML  method  of  [27],  for  hop 
timing  estimation  in  a  single-user  two-hop  setup  (one  hop  timing  instant  is  to  be  estimated),  which 
is  the  signal  model  considered  in  [27].  Then,  we  evaluate  the  performance  of  the  proposed  EM 
algorithm  in  more  general  cases. 


Performance  Comparison  for  Single-user  Two-hop  Scenario 

For  convenience  let  us  denote  the  iterative  ML  algorithm  of  [27]  as  IML.  We  compare  our  EM 
algorithm  to  IML.  Suppose  that  the  data  sequence  length  is  IV  =  128  and  the  signal  of  a  single 
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Figure  2:  Probability  of  detection  in  the  single-user  scenario. 

user  hops  once,  from  one  frequency  to  another,  in  this  data  sequence.  The  EM  algorithm  uses  a 
receive  array  with  M  —  6  antennas,  with  a  baseline  separation  of  half  a  wavelength  of  the  carrier 
frequency.  IML  uses  one  antenna.  Monte  Carlo  simulations  are  carried  out  to  compare  the  average 
performance  of  these  two  algorithms.  In  each  of  1000  realizations,  we  randomly  generate  not  only 
the  DOA,  hop  frequencies  (ranging  from  1  KHz  to  58  KHz)  and  complex  amplitudes,  but  also 
the  hop  instant  (ranging  from  the  30-th  snapshot  to  the  98-th  snapshot).  Note  that  in  [27]  the  hop 
instant  is  fixed  around  the  middle  of  the  data  sequence.  The  probabilities  of  detection  for  both 
algorithms  are  shown  in  Fig.  2  (b).  Here  if  the  estimated  hop  instant  is  the  same  as  the  true  hop 
instant,  the  detection  is  considered  successful.  From  Fig.  2  (b),  we  find  that  the  EM  algorithm 
significantly  outperforms  the  IMF  algorithm  at  low  SNR  range.  This  is  due  to  the  following  two 
reasons:  (i)  the  EM  algorithm  utilizes  multiple  antennas  while  IMF  uses  only  one  antenna;  (ii)  it 
is  not  shown  in  [27]  whether  IMF  can  guarantee  identifiability  for  frequency  estimation,  while  the 
MDF  algorithm  we  used  here  has  the  identifiability  guarantee. 

The  complexity  order  of  IMF  to  estimate  one  hop  instant  for  a  single  user  is  0(N2)  as  shown 
in  [27].  The  proposed  EM  algorithm  has  a  complexity  order  of  0(KMN 2)  for  multiple  signals 
with  K  hops.  Since  M  <§C  N,  we  conclude  that  the  complexities  of  the  two  approaches  for  one 
hop  instant  estimation  in  single  user  case  are  of  the  same  order. 

Multiple  Signals  with  Multiple  Hops 

We  evaluate  the  performance  of  the  EM  algorithm  for  the  case  of  multiple  signals  with  multiple 
hops.  The  receiver  array  has  M  =  6  antennas.  The  receiver’s  bandwidth  is  60  MHz.  Three  FH 
emitters  impinge  the  receiver  array  with  randomly  generated  DOAs  and  complex  amplitudes.  The 
hop  bandwidth  of  each  FH  emitter  is  80  MHz,  which  is  occupied  by  80  distinct  frequency  channels 
with  1  MHz  channel  spacing.  Therefore  there  is  a  transmitter-receiver  bandwidth  mismatch.  The 
data  size  is  6  x  400.  In  each  realization,  hop  instants  of  the  three  FH  signals  are  randomly  generated, 
and  hop  frequencies  are  randomly  chosen  from  the  80  channels. 

For  the  purpose  of  performance  evaluation,  a  successful  detection  of  hop  instant  is  defined  as 


8 


Figure  3:  The  EM  algorithm  for  hop  timing  estimation:  (a)  probability  of  detection;  (b)  probability 
of  false  alarm. 

follows:  for  a  given  true  hop  instant,  if  there  exists  an  estimated  hop  instant,  whose  deviation  from 
the  true  one  is  less  than  or  equal  to  4  snapshots,  then  this  detection  is  considered  successful.  For 
a  given  estimated  hop  instant,  if  its  distances  from  all  true  instants  are  greater  than  4  snapshots,  a 
false  alarm  is  claimed  to  have  occurred.  The  average  length  of  a  hop-free  segment  is  80  snapshots, 
therefore  4  snapshots  represent  5%  of  that  value.  Furthermore,  the  minimum  deviation  from  all 
estimated  hop  instants  to  a  given  true  hop  instant  is  defined  as  the  estimation  error  for  this  hop 
instant.  In  order  to  normalize  the  error,  we  divide  it  by  the  average  length  of  a  hop-free  segment. 

Fig.  3  (a)  plots  the  probability  of  detection  of  hop  instants  using  the  EM  algorithm,  where  both 
random  (“Rd-ini”)  and  spectrogram-based  (“Sp-Ini”)  initializations  are  tested.  The  corresponding 
probability  of  false  alarm  is  shown  in  Fig.  3  (b).  The  results  show  that  the  EM  algorithm  with 
spectrogram-based  initialization  does  a  good  job  for  hop  timing  estimation,  given  the  fact  that 
the  hop  sequences  and  model  order  variations  are  unknown.  It  can  be  seen  that  the  probability 
of  detection  is  about  95%  when  SNR  is  at  5  dB,  and  almost  100%  when  SNR  is  greater  than  15 
dB.  Correspondingly,  the  probability  of  false  alarm  is  less  than  5%  at  SNR=5  dB,  and  close  to 
zero  when  SNR  is  greater  than  15  dB.  The  EM  algorithm  with  spectrogram-based  initialization 
significantly  outperforms  its  randomly  initialized  counterpart. 

Because  of  using  a  simplified  EM  algorithm,  in  the  p-th  iteration,  the  k- th  hop  instant  n k  is 
only  searched  in  the  time  segment  bounded  by  n[^1  and  njj'  |l  i.  Hence  it  is  desirable  that  the 
initialization  is  approximately  in-line  with  each  system-wide  dwell,  which  can  be  achieved  by 
the  spectrogram  method  as  long  as  the  separation  of  two  adjacent  hop  instants  is  larger  than  the 
window  size  of  the  spectrogram. 

B.4.5  An  Application  Testbed  for  Bluetooth 

A  software  testbed  is  developed  to  apply  the  collision  resolution  method  in  multiple  co-located 
Bluetooth  piconets.  The  system  block  diagram  is  shown  in  Fig.  4.  Transmissions  from  multiple 
piconets  are  received  at  multiple  devices.  Without  loss  of  generality,  we  assume  that  each  piconet 
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Figure  4:  Collision  resolution  in  multiple  Bluetooth  piconets. 

consists  of  a  slave  and  a  master  device.  As  shown  in  Fig.  4,  a  user  in  piconet  1  may  receive  trans¬ 
missions  from  d  piconets.  The  transmissions  are  not  necessarily  synchronized  across  different 
piconets.  The  channel  is  assumed  to  be  flat-fading,  though  the  algorithm  can  be  extended  to  deal 
with  frequency-selective  channels.  The  proposed  receiver  consists  of  a  conventional  Bluetooth  re¬ 
ceiver  and  a  collision  resolution  unit.  We  assume  that  packet  collision  can  be  detected,  and  collided 
data  packets  are  sent  to  the  collision  resolution  unit  when  collisions  occur,  while  successful  packets 
are  sent  directly  to  the  output.  For  the  purpose  of  performance  comparison,  we  measure  the  BERs 
and  PERs  of  both  collision-resolution  enhanced  receiver  (BERe  and  PERe),  and  conventional  re¬ 
ceiver  (BERC  and  PERC).  The  software  testbed  for  the  collision  resolution  technique  is  developed 
based  on  Matlab  Simulink  environment  with  a  customized  graphic  user  interface.  It  simulates  data 
packet  transmissions  in  multiple  Bluetooth  piconets,  and  calculates  the  PER  and  BER.  A  number 
of  parameters  can  be  changed  in  the  testbed  to  simulate  difference  scenarios,  such  as  the  number 
of  piconets,  data  packet  length  and  SNR. 

B.5  Low  Complexity  Model  Order  Variation  Detection 

It  is  shown  in  the  previous  section  that  when  there  is  transmitter-receiver  bandwidth  mismatch, 
model  order  variation  has  to  be  detected  before  parametric  methods  can  be  applied.  Relying  on  a 
receiver  equipped  with  a  uniform  linear  array  (ULA),  we  have  developed  a  low-complexity  method 
for  the  detection  of  model  order  variations  in  the  context  of  FH  networks  with  transmitter-receiver 
bandwidth  mismatch  [41].  We  note  that  the  proposed  method  is  also  applicable  in  other  similar 
scenarios  where  the  transmitted  signals  are  uncorrelated. 

For  a  given  set  of  observations,  one  may  use  a  sliding-window  to  obtain  data  subsets,  and 
subsequently  apply  source  enumeration  techniques  such  as  SVD  or  MDL  to  each  subset  to  estimate 
its  model  order.  Model  order  change  points  can  then  be  detected  by  comparing  the  results  as  the 
data  window  slides.  However,  there  are  several  disadvantages  to  this  method.  First,  the  complexity 
of  SVD  and  MDL  is  high.  Second,  the  resolution  of  estimated  timing  is  limited  to  the  order  of  the 
window  size  at  low  SNR,  since  decisions  based  on  SVD  or  MDL  provide  hard-information  (the 
model  order).  Third,  SVD  and  MDL  break  down  when  the  number  of  signals  is  more  than  the 
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number  of  antennas. 

Here  we  also  adopt  a  sliding-window  approach.  More  importantly,  for  each  data  subset,  we 
compute  the  trace  of  the  signal  correlation  matrix  by  exploiting  inherent  data  structure.  The  trace 
sequence  provide  soft-information  that  later  is  fed  into  a  recursive  or  iterative  search  procedure 
to  obtain  accurate  timing  estimates  of  model  order  variations.  The  proposed  method  is  computa¬ 
tionally  much  simpler  than  sliding-window  SVD  or  MDL,  and  it  also  improves  the  probability  of 
detection  at  low  SNR.  Moreover,  the  method  works  even  if  the  number  of  effective  signals  is  more 
than  the  number  of  sensors. 

B.5.1  A  Trace  Calculation  Approach 

In  this  section  we  illustrate  how  soft-information  is  obtained  by  calculating  the  trace  of  the  corre¬ 
lation  matrix  of  each  data  subset.  Section  B.5.2  discusses  the  subsequent  detection  of  model  order 
variations  based  on  the  soft-information  sequence. 

The  data  model  is  given  in  (2).  Let  the  window  size  be  P.  We  assume  that  the  window  is 
small  enough  so  that  the  model  order  changes  at  most  once  in  each  window,  which  is  a  realistic 
assumption  for  slow  FH  systems,  and  if  the  sampling  rate  is  fast  enough  or  the  analysis  window  is 
short  enough,  this  assumption  is  also  valid  for  fast  FH  systems.  Without  loss  of  generality,  we  take 
a  data  subset 

X  [ajg  •  •  •  Xqpp_ ]J  , 

to  illustrate  the  trace  calculation  approach  in  the  following  two  cases:  constant  model  order  and 
varying  model  order  within  this  data  subset. 

Case  1:  Constant  Model  Order 

If  the  model  order  is  constant  and  equal  to  d  during  the  given  window,  the  M  x  P  data  matrix  can 
be  written  as 

X  =  AS  +  W, 

where  A  =  [/3ia(0i)  (ha{02)  ■  ■  ■  /3da(9d)\,  S  =  [s,  s,+i  •  •  •  sq+P_ i],  and  for  q  <  n  <  q  +  P, 
sn  =  \eJUJl’n'n  e3UJ2'n'n  ■  ■  ■  e]UJ,Ln '"]  .  W  is  the  corresponding  noise  matrix.  Suppose  the  noise  and 
signal  are  independent,  the  correlation  matrix  of  X  is  given  by 

Rx  =  E[XXh]  =  ARsAh  +  Rw,  (7) 

where  Rs  =  P\SSir]  is  the  signal  correlation  matrix,  and  the  noise  correlation  matrix  is  R\\  = 
a^I m-  Im  is  an  M  x  M  identity  matrix.  Assuming  that  signals  from  different  users  are  uncorre¬ 
lated,  the  signal  correlation  matrix  becomes  Rs  =  Id.  Hence  the  trace  of  Rx  is  given  by 

d 

tr (Rx)  =  tr (AhA)  +  Ma2w  =  M  ]T  |A|2  +  Ma2w ,  (8) 

i=  1 

where  we  have  used  the  fact  that  A  is  a  column-scaled  Vandermonde  matrix.  We  refer  to  the  trace 
of  Rx  as  the  “soft-information”,  which  is  a  linear  function  of  the  channel  power  as  shown  in  (8). 
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Notice  that  the  objective  here  is  to  detect  the  change  points,  rather  than  the  model  order  itself.  If  the 
change  points  are  known,  various  source  enumeration  techniques  can  be  applied,  e.g.,  SVD  [29], 
MDL  [31,32],  and  PDL  [40].  With  finite  samples,  Rs  ~  I <h  R-w  ~  cr^JM,  and  the  correlation 
matrix  of  X  can  be  estimated  as  Rx  =  p  YUn  ’q  '  xnxn- 

Case  2:  Varying  Model  Order 

We  now  illustrate  the  soft-information  in  the  presence  of  a  model  order  change  within  the  given 
data  subset  X.  For  simplicity,  assume  that  a  total  of  d  users  are  active  in  the  FH  system,  and  at 
time  q  +  l,  the  first  user  hops  out  of  the  observation  band.  Notice  that  if  more  than  one  users  hop 
out  at  the  same  time,  the  method  detects  the  model  order  variation  more  efficiently  because  of  the 
obviously  larger  change  in  (8).  We  can  write  the  data  matrix  as  X  =  AS  +  W ,  where 

&  [  Sq  '  Sq+l— 1  S q+l  '  '  '  &q+P—  1  ]  >  (9) 

and  for  n  <  q  +  l,  sn  is  the  same  as  in  Case  1,  and  for  n  >  q  +  l,  sn  is  the  same  as  in  Case  1  except 
that  its  first  element  is  zero.  In  this  case,  Rs  ^  Id.  In  fact,  it  is  shown  in  [41]  that  the  sample 
signal  correlation  matrix  is 

1 

Rs  =  p  SnSn  « 

n=q 

Similar  to  (8),  we  obtain 

tr(-Rx)  =  tr (ARsAh)  +  Ma2w  =  ^  7,  2  +  (H) 

i= 2 

It  can  be  seen  that  the  trace  of  the  estimated  correlation  matrix  with  base  at  n  =  q  is  linear  in  /, 
the  delay  or  lag  to  the  change  point.  If  there  are  multiple  changes  within  the  window,  the  situation 
is  more  complicated.  For  example,  an  incoming  signal  and  a  vanishing  signal  in  the  same  window 
could  cancel  each  other’s  effect  on  the  model  order.  Therefore  we  assume  that  the  window  size  is 
small  enough  that  the  model  order  changes  at  most  once  within  each  window. 

The  previous  derivation  assumes  that  A  is  a  Vandermonde  matrix.  If  A  is  not  a  Vandermonde 
matrix,  then  in  case  of  fixed  model  order,  (8)  becomes  tr  (Rx)  =  Ylt=  i  ||Oj||l  +  where  a*  is 
the  i-th  column  of  A,  and  1 1  •  ||2  stands  for  vector  2-norm.  In  the  case  of  varying  model  order,  (1 1) 
becomes  tr  (Rx)  =  ^||ai|||  +  ll^lli  +  ^aw-  Hence  the  trace  is  also  linear  to  the  “delay” 
to  the  change  point,  and  the  proposed  method  is  still  applicable. 

In  summary,  given  the  observations  x(n),  0  <  n  <  N,  we  can  obtain  a  soft-information 
sequence  by  using  a  sliding-window  as  follows  yq  =  tr  ^  xpxp^j ,  q  —  0, . . . ,  Q,  where  P 

is  the  window  size  and  Q  =  N  —  P.  Let  us  use  an  example  to  illustrate  y  =  [y0, . . . ,  yq]7  .  Suppose 
a  data  set  with  600  snapshots  is  available,  and  the  number  of  effective  signals  in  the  observation 
band  first  changes  from  three  to  one,  then  from  one  to  two,  due  to  frequency  hopping.  The  change 
points  are  170  and  410.  The  sequence  y  is  shown  in  Fig.  5,  where  the  window  size  is  50  and 


-  0 

p  w 


0  I 


d- 1 


(10) 


12 


7 


Figure  5:  Example  of  the  soft- information  sequence  obtained  by  trace  calculation  using  a  sliding- 
window,  where  true  points  of  model  order  changes  are  170  and  410. 

SNR  is  10  dB.  SNR  is  defined  as  10  log10(l/cr^).  In  Fig.  5,  we  observe  that  a  model  order  change 
manifests  itself  as  a  linear  variation  of  the  soft-information.  The  change  instants  are  given  by  the 
points  at  which  the  soft-information  flattens  out. 

Comparing  to  source  enumeration  algorithms  such  as  SVD  and  MDF,  the  trace  calculation 
method  is  of  much  lower  complexity,  and  works  even  when  the  number  of  the  signals  is  more  than 
the  number  of  antennas.  For  example,  for  an  M  x  P  data  subset,  the  complexity  is  about  O(MP) 
for  trace  calculation,  and  0(M2P  +  P3)  for  SVD  and  MDF.  After  obtaining  the  soft-information 
sequence  y ,  the  next  step  is  to  estimate  the  change  points,  rik,  for  k  =  0,  •  •  •  ,  K  —  1.  This  is 
discussed  next. 

B.5.2  Change  Detection  Methods 

We  have  developed  three  methods  to  estimate  the  change  points  from  a  sequence  y ,  which  can 
be  a  soft-information  sequence  obtained  by  the  trace  calculation  or  a  hard-information  sequence 
obtained  by  SVD  or  MDF.  One  method  is  a  recursive  search  similar  to  the  dynamic  programming 
principle  (see,  e.g.,  [42,  Chapter  12]),  and  the  other  is  a  decoupled  iterative  search  similar  to  the 
expectation  maximization  principle  (see,  e.g.,  [43]).  However,  we  note  that  since  the  detection  is 
not  maximum  likelihood,  and  we  do  not  pursue  dynamic  programming  nor  expectation  maximiza¬ 
tion  algorithms.  The  third  method  is  a  simple  low-complexity  difference  approach  to  estimate  the 
change  points  from  the  soft-information  sequence.  More  details  can  be  found  in  [41]. 

B.5.3  Simulation  Results 

We  present  simulation  results  to  illustrate  the  performance  of  the  proposed  methods.  For  a  given 
data  matrix  with  model  order  variations,  we  apply  the  trace  calculation  (TC),  SVD,  and  MDF  ap¬ 
proaches  to  data  subsets  obtained  using  a  sliding-window.  Subsequently  the  recursive,  the  iterative 
and  the  difference  detection  methods  are  used  to  estimate  the  change  points.  Notice  that  SVD  and 
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Table  1:  Probability  of  detection  of  model  order  variation  (AWGN  channels) 


SNR 

OdB 

3  dB 

6  dB 

9  dB 

Pd 

% 

SVD-R 

3 

28 

89 

98 

SVD-I 

3 

28 

87 

98 

MDL-R 

7 

42 

91 

99 

MDL-I 

6 

42 

90 

99 

TC-R 

51 

71 

91 

97 

TC-I 

36 

53 

81 

94 

TC-D 

82 

91 

97 

98 

Table  2:  Probability  of  detection  of  model  order  variation  (flat-fading  channels) 


SNR 

OdB 

3  dB 

6  dB 

9  dB 

Pd 

% 

SVD-R 

33 

51 

83 

95 

SVD-I 

30 

50 

84 

94 

MDL-R 

35 

55 

87 

95 

MDL-I 

33 

55 

86 

95 

TC-R 

59 

72 

84 

91 

TC-I 

46 

59 

75 

85 

TC-D 

80 

87 

92 

94 

MDL  return  hard-information  -  the  actual  decision  on  model  orders  for  each  data  subset.  The 
decision  threshold  used  in  the  SVD  method  is  chosen  according  to  the  modified  third  bound  given 
in  [29].  The  MDL  algorithm  used  here  follows  [31]. 

Monte  Carlo  simulations  are  conducted  to  assess  the  performance  of  the  proposed  method. 
The  sources  are  three  narrow-band  FH  emitters  with  randomly  generated  DOAs.  We  generate  200 
realizations  and  compute  the  probability  of  detection  for  each  SNR.  In  each  realization,  the  model 
order  change  points  and  the  number  of  signals  that  hop  in/out  of  the  observation  band  are  randomly 
generated.  A  detection  is  defined  successful  when  the  difference  between  an  estimated  instant  and 
the  true  one  is  less  than  or  equal  to  5,  which  is  ten  percent  of  the  window  length  50.  Table  1 
shows  the  probability  of  detection  ( Pd )  using  the  recursive  method  (R),  the  iterative  method  (I) 
and  the  difference  method  (D)  for  AWGN  channels,  where  all  /3^’s  are  fixed  to  1.  Table  2  shows 
the  probability  of  detection  using  these  three  methods  for  flat-fading  channels,  where  all  f5^k  s  are 
randomly  generated  with  standard  complex  Gaussian  distribution. 

For  both  recursive  and  iterative  detections,  the  results  indicate  that  the  TC  based  approach 
outperforms  the  SVD  and  MDL  based  methods  at  low  SNR  regime.  This  is  mainly  because  that 
at  low  SNR,  the  resolution  of  SVD  and  MDL  is  limited  by  the  window  size,  while  TC  does  not 
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have  this  problem  since  it  returns  soft-information.  The  iterative  detector  slightly  underperforms 
the  recursive  detector  but  comes  with  much  lower  complexity.  The  difference  method  combined 
with  trace  calculation  has  the  least  complexity  and  competitive  performance.  When  SNR  is  greater 
than  10  dB,  all  three  methods  (SVD,  MDL,  and  TC)  can  achieve  a  probability  of  detection  close 
or  equal  to  100%. 


B.6  Optimizing  Eigenvector-Based  Frequency  Estimation 

Multidimensional  frequency  estimation  plays  an  important  role  in  our  proposed  interference  miti¬ 
gation  framework.  Recently  an  eigenvector-based  algorithm  has  been  developed  for  multidimen¬ 
sional  frequency  estimation  with  a  single  snapshot  of  data  mixture.  Unlike  most  existing  algebraic 
approaches  that  estimate  frequencies  from  eigenvalues,  the  eigenvector-based  algorithm  achieves 
automatic  frequency  pairing  without  joint  diagonalization  of  multiple  matrices,  but  it  fails  when 
there  exist  identical  frequencies  in  certain  dimensions  because  eigenvectors  are  not  linearly  in¬ 
dependent  anymore.  We  develop  an  eigenvector-based  algorithm  for  multidimensional  frequency 
estimation  with  finite  data  snapshots.  We  introduce  complex  weighting  factors  so  that  the  algorithm 
is  still  operational  when  there  exist  identical  frequencies  in  one  or  more  dimensions.  Furthermore, 
the  weighting  factors  are  optimized  to  minimize  the  mean  square  errors  of  the  frequency  estimates. 
Simulation  results  show  that  the  proposed  algorithm  offers  competitive  performance  when  com¬ 
pared  to  existing  algebraic  approaches  but  at  lower  complexity. 

The  data  model  for  iV-D  frequency  estimation  can  take  a  variety  of  forms.  Here  we  refer  to  a 
single  snapshot  iV-D  frequency  mixture  as  an  iV-D  array  X  with  typical  element 


•Uni, m2,--  ,mjv 


F  N 


/= 1  n=  1 


1) 


+  W 


mi 


(12) 


where  mn  =  1, . . . ,  Mn,  for  n  —  1, . . . ,  N,  and  Mn  is  the  dimension  size  of  the  n-th  dimension. 
The  total  sample  size  is  M  :=  Y\n=i  Mn-  In  (12),  the  frequencies  tOfyn  G  (— tt,  i r],  for  /  =  1, . . . ,  F, 
n  —  1, . . . ,  N,  and  wm  1,m2,-,mN  is  white  Gaussian  noise  with  variance  cr'2.  Similarly  T  snapshots 
of  iV-D  frequency  data  mixtures  may  be  modeled  as  T  N- D  arrays,  X  (t  ).  with  typical  element 


F  N 

xmi,m2,...  ,mN(t)  =  y ^Cf(t)  TT  +wmi|ma,...,mAf(t),  t  =  l,--  -  ,  T,  (13) 

/= 1  n= 1 

where  t  is  the  snapshot  index,  which  can  be  a  time  index,  or  trial  index  in  case  of  multiple  trials 
of  experiments.  T  =  1  corresponds  to  the  single  snapshot  case.  The  frequency  set  (c <Jf,n}n=i  1S 
aiV-D  frequency  component,  and  there  are  F  such  components.  The  objective  of  N- D  frequency 
estimation  is  to  estimate  {uf>n}^=1,  for  /  =  1, . . . ,  F,  from  given  X(t),  t  —  1,  •  •  •  ,  T.  Notice  that 
the  same  data  model  for  multiple  snapshot  case  has  been  used  in  [44,45]. 

In  order  to  facilitate  the  presentation,  we  introduce  the  equivalent  data  models  based  on  Khatri- 
Rao  product.  Given  (13),  define  the  sample  vector  x(t)  for  the  t-th  snapshot  as 

x{t)  =  ./•|.i.....i(/)  .r,.i 2(t)  •••  Xi.i,. ..,2,1  {t)  •••  XMl ,M2,...,MJV(f)]T 
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Furthermore,  define  N  Vandermonde  matrices  An  G  cMnXF  with  generators  = ,  such  that 

An  :=  [ahn  a2,n  ■  ■  ■  aF,n] ,  where  a/>n  =  [l  eJU}f’n  ■  ■  ■  ej(Mn~1)u}f’n] T  ,  n  =  1, . . . ,  N. 

It  can  be  verified  that 

x(t)  —  Ac(t)  +  w(t),  t  =  l,...,T,  (14) 

where  w(t)  is  the  noise  vector,  and 

A  :=  Ai  ©  A-2  ©  •  •  •  ©  An ,  c{t)  :=  [ci(t)  c2(t )  •  •  •  cF(t)]T . 

Define 

X  :=  [*(1)  a;(2)  •  •  •  *(T)]  G  CMxT,  C  :=  [c(l)  c(2)  •  •  •  c(T )]  G  CFxT, 

then  the  data  model  in  (14)  can  be  rewritten  in  matrix  form  as 

X  =  AC  +  W,  (15) 

where  W  is  the  corresponding  noise  matrix.  We  will  need  the  following  lemmas  [46]. 


Lemma  1  Define  a  set  of  N-D  selection  matrices  as 


:=  Jf1  ®  J*2 


Jh/2,-,eN  ■=  Jf1  ®  Ji*  ■  ■  ■  ®  Ji^,  (16) 

J<Ln  :=  [°^x(4-i)  Iku  0Knx(Ln-en)\  ,  (17) 

where  ln  =  1, . . . ,  Ln,  and  Kn  and  Ln  are  positive  integers  satisfying  Kn  +  Ln  =  Mn  +  1  for 
n  —  1 , ,N.  Further  define  an  N-D  smoothing  operator  for  the  snapshot  vector  in  (14)  as 

S[x(t)}  ■=  [Ji,!,...  ,!*(<)  j  1,1,...  px(t)  J  ,LNx(t)  ix(t)  ■■■  JLl,L2,-,LNx(t)]  , 

then  it  can  be  verified  that  in  the  absence  of  noise 

Xs{t)  :=  S[x(t)]  =  (Afl]  ©  Af2)  ©  •  •  •  ©  A^})L>(c(f))  (aSLi)  ©  A[L2)  ©  •  •  •  © 

Lemma  2  Given  N  Vandermonde  matrices  An  G  CM" x F,  with  generators  for  n  = 

1, . . . ,  N,  and  a  complex  matrix  C  G  CFxT,  if  we  define 

B  :=  _  n tchd((3)  _  ’  (18) 

where  IIT  is  aT  x  T  permutation  matrix  with  ones  on  its  anti-diagonal,  and 


f3  =  [e~j^  e~jf)2  ■  ■  ■ 

N 

fif  =  y,  (Mn  -  1  V/,n, 


~30f]T 


then  the  rank  of  the  matrix 

H  :=  A[Ll)  ©  Ail2)  •  •  •  ©  A^n)  ©  B  (20) 

is  min  |2 T  Hn=i  7-n,  F  |  almost  surely,  provided  that  the  N F  frequencies  {wf,n}n=i>  /  =  1,  •  •  •  ■  F, 
and  the  TF  elements  of  C,  are  drawn  from  distributions  that  are  continuous  with  respect  to  the 
Lebesgue  measure  in  QNFxl  and  CTFxl,  respectively. 
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B.6.1  The  Eigenvector-Based  Algorithm  for  Ar-D  Frequency  Estimation 

In  this  section  we  present  the  algorithm  for  iV-D  frequency  estimation  from  multiple  snapshots. 
For  simplicity  of  exposition,  the  algorithm  is  developed  in  the  noiseless  case.  The  noise  effect  on 
the  performance  of  the  algorithm  is  analyzed  in  Section  B.6.2. 

Given  (14)  in  the  noiseless  case,  we  can  apply  the  smooth  operator  S  defined  in  Lemma  1  to 
every  snapshot  x(t),  and  obtain 

Xs(t)  ■=  S[*(t)]  =  GD(c(t))HT,  (21) 

where 

G  :=  A[Kl)  ©  Af2)  •  •  •  0  A{*N\  H  :=  A[Ll)  0  Af- 2)  •  •  •  ©  A^N\ 

The  positive  integers  Kn  and  Ln,  n  —  1,  •  •  •  ,  N,  are  chosen  such  that 

Kn  +  Ln  =  Mn  +  1,  1  <  n  <  N.  (22) 

To  further  explore  the  data  structure,  we  can  perform  the  forward-backward  smoothing  on  the  data 
vector  x(t)  in  (14).  Define  y(t)  :=  II MX*(t),  where  II m  is  an  MxM  permutation  matrix  with 
ones  on  its  anti-diagonal.  It  can  be  verified  that?/  (t)  =  Ac(t),  where  c(t)  =  \ci(t),  c2(f),  •  •  •  ,Ci?(t)]T, 
with  Cf(t )  =  ,  and  /3f  is  defined  in  (19).  Applying  the  same  technique  to  y(t)  that  we 

used  to  construct  Xs(t)  from  x(t),  we  obtain 

Ys(t)  :=S[y(t)]=GD(c(t))HT. 

We  then  collect  all  the  smoothed  data  matrices  to  obtain 

X:=[XS{1)  Xs{2)  ■  ■■  XS{T)  Y S(T )  Y S(T  -  1)  ^(1)].  (23) 

It  can  be  verified  that 

X  =  G(H  ©  B)t  =  GHT,  (24) 

where  B  and  H  are  defined  in  (18)  and  (20),  respectively.  A  key  step  of  our  algorithm  is  the  con¬ 
struction  of  X  to  ensure  that  it  is  of  rank  F  almost  surely.  Note  that  similar  smoothing  technique 
has  been  used  in  [44],  but  its  relation  to  the  Khatri-Rao  product  is  not  explored.  In  (24),  since  G 
is  the  Khatri-Rao  product  of  multiple  Vandermonde  matrices,  G  is  almost  surely  full  column  rank 
if  nli  Kn  >  F.  According  to  Lemma  2,  if  2 T  nli  Ln  >  F,  H  has  full  column  rank  almost 
surely.  According  to  the  Sylvester’s  inequality  [47] 

- — X1  - — -X  ' — X 

rank(G)  -frank  (IT  )  —  F  <  rank  (GIT  )  <  min{rank(G), rank(TT  )}, 
hence  X  is  of  rank  F  almost  surely.  The  singular  value  decomposition  (SVD)  of  X  yields 

X  =  USYSV^,  (25) 
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where  Us  has  F  columns  that  together  span  the  column  space  of  X.  Since  the  same  space  is 
spanned  by  the  columns  of  G ,  there  exists  an  F  x  F  nonsingular  matrix  T  1  such  that 

Us  =  GT  1  (26) 

Similar  to  the  IMDF  algorithm  [48],  once  the  signal  subspace  Us  is  obtained,  we  can  construct 
two  matrices  from  U s  whose  general  eigenvalues  are  the  exponentials  of  the  first  dimension,  and 
then  use  the  general  eigenvectors  to  estimate  N-D  frequencies.  However,  as  mentioned  before, 
the  IMDF  algorithm  fails  when  there  exist  identical  frequencies  in  the  first  dimension  since  the 
eigenvectors  are  not  linearly  independent  anymore.  Furthermore,  it  has  been  shown  in  [48]  that 
the  performance  of  the  IMDF  algorithm  is  severely  degraded  if  there  are  close  frequencies  in 
the  first  dimension.  To  address  this  problem,  in  the  following,  we  present  a  method  to  construct 
two  matrices  whose  general  eigenvalues  are  weighted  sum  of  the  N-D  exponentials.  The  N- D 
frequencies  are  still  resolved  from  the  general  eigenvectors. 

We  define  two  selection  matrices  J  \  and  J2  as 


TV 

J 1  J  1,\  ®  J  1,2  '■  '  ®  J  1,7V,  J 2  =  OinJ 2, re,  (27) 

n= 1 

where  J \,n  =  0(Kn-i)xi],  J2,n  =  [0(ic„-i)xi  Ik„-i]  ,  J 2,n  =  (•T'i.i  ®  •••<£>  «7i,n-i)  ® 

J 2 ,n  ©  (J l,n+l  ®  ®  J l,7v)- 

Here  {an}%=1  are  complex  weighting  factors,  which  can  be  randomly  chosen  initially.  As  we 
will  show  in  Section  B.6.2,  the  MSEs  of  the  frequency  estimates  are  affected  by  these  weighting 
factors  in  the  noisy  case.  Next,  we  obtain  two  equal-sized  matrices  U  \  and  U2  by 

Ul  :=  J,US,  U2  :=  JoUs.  (28) 

According  to  the  property  of  Khatri-Rao  product  [49]:  ( A  0  B)(C  G  D)  =  AC  ©  BD ,  it  can 
be  verified  that 

Ui  =  PT  \  U2  =  PD(C)T-\  (29) 

where  P  =  A\K'^V)  ©  At2R'2^1'1  0  •  •  •  0  It  is  clear  that  P  has  full  column  rank  almost 

surely  if  F  <  f\^=1(Kn  -  1).  In  (29),  C  :=  [Ci,  C2,  •  •  •  ,  C f}t,  and  (f  =  J2n=i  aneJUJf’n.  Therefore 
we  have 

U\U  2  =  TD{C)T-\ 

Clearly  we  can  choose  {an}^=1  to  ensure  the  elements  of  C  are  distinct  even  if  there  exist  identical 
frequencies  in  one  or  more  dimensions,  but  this  is  not  guaranteed  by  randomly  generated  {an}^=1. 
We  will  discuss  how  to  choose  the  weighting  factors  in  Section  B.6.2.  T  can  be  retrieved  from 
the  eigen-decomposition  of  U\U2  up  to  column  permutation  and  scaling  ambiguity.  Suppose  that 
the  eigen-decomposition  of  U\U2  gives  Tsp  =  TAA,  where  A  is  a  nonsingular  diagonal  column 
scaling  matrix  and  A  is  a  permutation  matrix.  Once  we  obtain  Tsp,  we  can  retrieve  P  up  to  column 
permutation  and  scaling  ambiguity  according  to 

Psp  =  U 1  Tsp  =  PAA.  (30) 
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Notice  that  P  is  the  Khatri-Rao  product  of  N  Vandermonde  matrices,  and  there  are  F  columns  in 
P.  The  N  frequencies  of  the  same  N- D  component  appear  in  the  same  column  of  P.  In  other 
words,  for  fixed  /,  (c ■Of,n}n=i  appear  in  the  same  column  of  P.  Thanks  to  this  structure,  we  can 
obtain  F  N- D  frequency  components  by  dividing  suitably  chosen  elements  of  the  aforementioned 
columns  of  Psp.  Therefore  the  column  scaling  and  permutation  will  not  have  a  material  effect 
on  the  algorithm.  For  this  reason,  we  may  drop  subscript  “sp”  from  now  on  as  long  as  it  is  clear 
from  the  context.  Suppose  appear  in  the  /- th  column  of  P,  then  each  of  them  can  be 

obtained  by  anyone  of  the  following  quotients 

e?“f'n  =  mod  (k  -  1,  <_,)  >  K'n,  for  /  =  1, . . . ,  F,  (31) 

Pk-K'n,f 


where  1  <  k  <  K'0,  pkj  is  the  (k.  /)- th  element  of  P,  and 


K 


n;L+1(Ap-l).  0<n<N  —  1, 

1,  n  =  N. 


(32) 


Notice  that  the  frequencies  are  automatically  paired  because  the  frequencies  n  —  1 , . . .  ,N) 
of  the  same  A’-D  component  (the  /-th  component)  are  obtained  from  the  same  column  of  P. 

If  the  data  observations  are  noisy  as  given  in  (14),  applying  the  above  algorithm  we  can  obtain 
the  estimate  of  P  as  P.  In  order  to  reduce  the  MSEs  of  frequency  estimates,  we  use  the  average 
of  all  the  quotients  in  (31)  to  obtain  an  estimate  of  the  Ar-D  exponential.  Therefore,  can  be 
estimated  by  the  following  average 

e*"/."  =  —  V  /kJ  ,  n  =  1, . . . ,  N,  (33) 

Fn  j~t  Pk-KU 

mod  (  fc  —  1 ,  Krn  _  j  )  >  K*n 


where  fin  =  K'0(Kn  —  2)/ (Kn  —  1).  The  average  is  also  the  so-called  “circular  mean”  in  direction 
statistics  [50].  Finally  the  frequency  estimates  are  obtained  by 

w/,„  =  X  (log  e**’/.")  .  (34) 

After  the  frequency  estimates  are  obtained,  the  amplitude  matrix  C  can  be  obtained  by  solving 
(15)  using  a  least-squares  approach. 

For  an  Ar-D  frequency  estimation  algorithm,  the  maximum  number  of  uniquely  resolvable 
frequency  components  in  the  absence  of  noise  is  referred  to  as  its  identi fi ability  bound.  We  sum¬ 
marize  our  main  result  on  statistical  identifiability  for  the  proposed  algorithm  in  the  following 
theorem  [46]. 


Theorem  1  Given  T  snapshots  of  sums  of  F  N-D  exponentials  as  in  (13),  in  the  absence  of  noise, 
the  parameter  set  {cf(t)}I=  i)>  f  —  1,  •  •  • ,  F,  is  almost  surely  uniquely  identifiable  by 

the  proposed  algorithm  in  Section  B.6.1,  provided  that 


N 


F  < 


max  min  1 

n 

Kn-\-Ln=Mn-\- 1 

1  <Kn<Mn 

\n.=  l 

l<n<N 
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N 


n=  1 


(35) 


where  the  NF  frequencies  {wf,n}n=i>  f  =  k  •  •  •  ?  7b  and  TF  complex  amplitudes  {cf(t)}J=1, 
f  =  1 , ,F,  are  assumed  to  be  drawn  from  distributions  that  are  continuous  with  respect  to  the 
Lebesgue  measure  in  0NFx  1  and  CTFxl,  respectively. 


B.6.2  Optimization  of  the  Eigenvector-Based  Algorithm 


We  have  derived  the  theoretic  error  variance  of  the  eigenvector-based  algorithm  in  [46] .  We  have 
also  derived  the  Cramer-Rao  Bound  for  multidimensional  frequency  estimation  models  [46].  Sup¬ 
pose  that  Uftn  —  u 0f,n  +  Ac u/>n.  It  is  shown  in  [46]  that 


V 


E[Aw2f  ] 

/  /  llln  - 7^ - 7 

ct2^°  varCRB 


(X  7(a), 


(36) 


where  7 (a)  =  Yf}=l  Ej=/+i  | and  U>f  =  Uf’2  ' ' '  Uf’N ^  T°  minimize  the  error 
variance,  an  optimal  a  can  be  obtained  by 


aopt  =  arg min 7(a),  subject  to  ||a||<l. 


(37) 


The  optimization  problem  (37)  is  a  so  called  sum-of-ratios  fractional  programming  problem, 
which  is  a  difficult  global  optimization  problem  [51].  There  is  no  efficient  algorithm  available  to 
solve  it  to  date.  We  propose  to  use  grid  search  in  the  super-sphere  ||a||  <1  to  find  a  moderate  initial 
value  of  a,  then  use  a  Newton  type  algorithm  to  find  an  optimal  {an}%=1.  In  order  to  reduce  the 

complexity,  we  may  set  \an\  —  for  n  =  1,  •  •  •  ,  N,  and  the  search  grid  does  not  need  to  be 
fine  (for  example,  the  step  size  of  angle  in  one  dimension  can  be  set  to  7 t/F). 

Alternatively  we  can  use  the  following  method  to  get  a  moderate  initial  value  of  {an}^=1.  If 
we  define 

(38) 


eial  :=  mm 


then  we  have 


7(a)  < 


F(F  —  1) 


e  a 


If  we  can  solve  the  following  optimization  problem 


a0  =  arg  max  e(a),  (39) 

l|a||<i 

the  upper  bound  of  7(a)  is  minimized.  The  optimization  problem  (39)  can  be  solved  using  a  se¬ 
quential  quadratic  programming  (SQP)  method,  which  is  a  common  quasi-Newton  type  algorithm 
available  in  many  optimization  packages  such  as  the  optimization  toolbox  in  Matlab. 

The  algorithm  using  an  optimal  a  for  A’-D  frequency  estimation  is  described  in  Table  3.  Notice 
that  the  first  step,  which  involves  the  SVD  of  X ,  is  most  computationally  complex.  But  this 
step  is  executed  only  once.  To  appreciate  the  proposed  algorithm  in  Table  3,  we  compare  three 
methods  to  obtain  an  optimal  a.  The  difference  is  only  in  Step  3  of  Table  3,  where  we  may:  (a) 


20 


Table  3:  An  improved  eigenvector-based  algorithm  using  optimal  weighting  factors 


1.  Given  (14),  follow  (21)— (25)  to  obtain  U s. 

2.  Randomly  select  a  subject  to  \an\  =  n  =  1, . . . ,  N, 
compute  {uJf,n}n=i^  f  =  1,  ■  •  • ,  F,  using  (27)-(34). 

3.  Based  on  {cD/,n}^=1>  for  /  —  1 .... .  /•',  obtain  an  updated  aopt 
by  first  solving  (39)  using  SQP  to  get  initials,  then  solving  (37) 
using  a  Newton  method. 

4.  Compute  updated  {w/,n}^=1,  1  <  /  <  F,  with  aopt  using 
(27M34). 

5.  Iterate  Steps  3  and  4  until  frequency  estimates  converge  (typ¬ 
ically  one  execution  of  Steps  3-4  is  sufficient). 


using  the  solution  of  (39)  as  aopt ,  referred  to  as  “Minmax”  here;  (b)  using  the  solution  of  (39)  as 
initials,  then  solving  (37)  to  obtain  raopt,  referred  to  as  “Minmax+Newton”,  which  is  the  proposed 
algorithm;  and  (c)  solving  (37)  by  grid  search  first  then  refine  it  using  a  Newton  method,  referred 
to  as  “Grid+Newton”.  These  approaches  are  applied  to  estimate  three  3-D  frequency  components 
from  3  snapshots  of  6  x  6  x  6  noisy  data  samples.  Fig.  6  (a)  depicts  the  Root  Mean  Square 
Error  (RMSE)  versus  Signal-to-Noise  Ratio  (SNR).  The  RMSE  is  obtained  by  averaging  over  all 
frequencies  after  only  one  iteration  of  Steps  3-4,  except  for  the  case  indicated  with  “Random  a” 
where  only  Steps  1-2  are  executed.  The  corresponding  CRB  on  STD  is  also  plotted. 

It  can  be  seen  from  Fig.  6  (a)  that  the  three  optimization  methods  provide  similar  perfor¬ 
mance,  and  all  outperform  the  case  with  randomly  chosen  «.  It  turns  out  one  iteration  of  Steps 
3-4  is  sufficient,  as  demonstrated  by  Fig.  6  (b),  where  we  plot  the  RMSE  of  frequency  esti¬ 
mates  versus  the  number  of  iterations  of  Steps  3-4.  Zero  iteration  corresponds  to  the  case  with 
only  randomly  chosen  ex.  We  observe  that  one  execution  of  Steps  3  and  4  is  sufficient  as  further 
iterations  only  provide  negligible  performance  improvement  if  any.  It  can  be  seen  that  “Min- 
max+Newton”  and  “Grid+Newton”  are  comparable,  and  both  are  slightly  better  than  “Minmax”. 
Because  “Grid+Newton”  has  a  higher  complexity  than  “Minmax+Newton”,  we  choose  “Min¬ 
max+Newton”  with  one  iteration  as  the  proposed  algorithm  in  Table  3  for  optimizing  the  weighting 
factors,  which  is  the  algorithm  used  in  the  simulations  of  Section  B.6.3. 

B.6.3  Simulation  Results 

In  this  section  we  present  the  Monte  Carlo  simulation  results  to  demonstrate  the  performance 
(measured  by  RMSE)  of  the  optimized  eigenvector-based  frequency  estimation  algorithm,  which 
is  also  compared  to  other  Ar-D  frequency  estimation  algorithms  as  well  as  the  associated  CRB. 
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Figure  6:  (a)  RMSE  of  different  optimization  methods  versus  SNR;  (b)  RMSE  of  different  opti¬ 
mization  methods  versus  the  number  of  iterations  of  Steps  3-4. 

2-D  Identical  Frequency  Estimation  from  Single  Snapshot 

In  the  first  experiment,  the  proposed  algorithm,  MEMP  [37],  Unitary  ESPRIT  [44]  and  2-D  ES¬ 
PRIT  [34]  are  applied  to  estimate  three  2-D  frequency  components  from  a  20x20  noisy  data  set. 
The  amplitudes,  c/(  1)  for  /  =  1, . . . ,  F,  are  set  to  be  one  for  this  case.  The  three  frequency  pairs 
are  cuii2)  =  (0.557T,  0.207t),  (cu2,i,  1^2,2)  =  (0.607T,  0.207t),  and  (u;3ji,  0*3,2)  =  (0.607T,  0.257t). 
Notice  that  there  are  identical  frequencies  in  both  dimensions,  which  is  a  case  that  the  IMDF  al¬ 
gorithm  fails  to  deal  with.  Fig.  7  depicts  the  performance  comparison.  In  Fig.  7  (a),  we  plot  the 
RMSE  of  various  algorithms  and  the  average  CRB  on  STD  in  the  two  dimensions.  The  RMSE 
results  are  averaged  over  all  frequencies  and  obtained  through  1000  realizations.  In  Fig.  7  (a), 
“Proposed  algorithm”  refers  to  the  one  with  one  iteration  of  Steps  3-4  using  “Minmax+Newton”. 
For  our  proposed  algorithm,  the  smoothing  parameters  ({Kn}^=1,  {Ln}2n=1)  are  chosen  such  that 
the  identifiability  bound  can  be  achieved.  The  parameters  for  other  algorithms  are  chosen  accord¬ 
ing  to  their  respective  references.  As  shown  in  Fig.  7  (a),  the  proposed  algorithm  offers  comparable 
performance  as  that  of  Unitary  ESPRIT,  and  outperforms  2-D  ESPRIT  and  MEMP. 

In  Fig.  7  (b),  we  compare  the  optimized  weighting  factors  to  randomly  chosen  ones,  where 
“Random  a ”  means  zero  iteration  of  Steps  3-4  in  Table  3.  The  theoretic  RMSE  is  obtained  with 
an  optimized  a.  by  solving  (37)  using  the  true  frequencies,  which  serves  as  a  benchmark  since  in 
our  algorithm  ol  is  optimized  when  the  true  frequencies  are  unknown.  It  is  clear  that  the  proposed 
algorithm  significantly  outperforms  the  one  with  random  weighting  factors,  and  the  simulated 
RMSE  of  the  proposed  algorithm  matches  well  to  the  theoretic  RMSE  for  moderate  to  high  SNR. 

2-D  Close  Frequency  Estimation  from  Multiple  Snapshots 

In  the  second  experiment,  the  proposed  algorithm.  Unitary  ESPRIT  and  RARE  are  applied  to 
estimate  three  2-D  frequency  components  from  10  snapshots  of  noisy  data,  each  of  size  12  x  12,  as 
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Figure  7:  (a)  Comparison  of  different  algorithms  for  2-D  frequency  estimation  from  single  snap¬ 
shot;  (b)  Comparison  of  optimized  at  and  randomly  chosen  at. 

given  in  (14).  The  amplitudes,  c/(f),  for  /  =  1, . . . ,  F  and  t  —  1, _ ,  T,  are  drawn  from  a  complex 

Gaussian  distribution.  The  three  frequency  pairs  are  (u^i,  =  (0.727T,  0.627t),  {u2, i ,  6^2,2)  = 

(0.747T,  0.587t),  and  (cu3ji,  cc3i2)  =  (0.767T,  0.607t).  Notice  that  frequencies  are  close  to  each  other 
in  both  dimensions.  Fig.  8  depicts  the  simulated  RMSE  of  various  algorithms,  along  with  the 
corresponding  CRB  and  the  theoretic  RMSE  of  the  proposed  algorithm.  Multidimensional  data 
smoothing  is  also  performed  for  the  Unitary  ESPRIT  algorithm  and  the  RARE  algorithm.  It  can  be 
seen  from  Fig.  8  (a)  that  the  proposed  algorithm  offers  competitive  performance  when  compared 
with  the  Unitary  ESPRIT  and  RARE  algorithms.  Note  that  the  proposed  algorithm  has  lower 
complexity  than  these  two  algorithms,  because  the  proposed  algorithm  does  not  require  a  frequency 
pairing  step  (the  Unitary  ESPRIT  algorithm  achieves  automatic  frequency  pairing  through  iterative 
joint  diagonalization).  Fig.  8  (b)  confirms  again  that  optimized  weighting  factors  outperform 
randomly  chosen  weighting  factors,  and  the  simulated  RMSE  of  the  proposed  algorithm  matches 
its  theoretic  RMSE  at  high  SNR. 

3-D  Identical  Frequency  Estimation  from  Multiple  Snapshots 

In  the  third  experiment,  the  proposed  algorithm  and  the  Unitary  ESPRIT  algorithm  are  applied  to 
estimate  three  3-D  frequency  components  from  10  snapshots  of  6x6x6  noisy  data  samples.  There 
are  identical  frequencies  in  all  dimensions.  The  amplitudes  are  drawn  from  a  complex  Gaussian 
distribution.  Fig.  9  shows  the  performance  comparisons.  From  Fig.  9,  we  notice  the  proposed 
algorithm  also  offer  competitive  performance  in  3-D  frequencies  estimation  compared  to  the  Uni¬ 
tary  ESPRIT  algorithm.  Because  the  pairing  strategy  of  RARE  algorithm  is  not  applicable  when  all 
three  dimensions  have  identical  frequencies,  we  do  not  include  RARE  in  this  experiment.  Notice 
that  the  simulated  RMSE  of  the  proposed  algorithm  matches  its  theoretic  RMSE  for  moderate  to 
high  SNR  range. 


23 


Figure  8:  (a)  Comparison  of  different  algorithms  for  2-D  frequency  estimation  from  multiple  snap¬ 
shots;  (b)  Comparison  of  optimized  «  and  randomly  chosen  «. 

B.7  Conclusions 

In  this  project  we  have  proposed  a  signal  processing  framework  for  co-channel  interference  mit¬ 
igation  when  multiple  FH  network  coexist.  The  EM  algorithm  jointly  estimates  hop  timing  and 
frequency  estimation  of  multiple  FH  signals,  with  unknown  hop  sequences  and  possible  bandwidth 
mismatch.  A  simple  initialization  step  based  on  the  data  spectrogram  gives  initial  estimates  of  hop 
timing  for  the  EM  algorithm.  Simulation  results  have  shown  that  the  EM  algorithm  is  capable  of 
obtaining  the  operation  characteristic  of  noncooperative  FH  emitters.  When  there  is  a  transmitter- 
receiver  bandwidth  mismatch  in  multiple  FH  networks,  the  model  order  changes  evenif  the  number 
of  active  emitters  does  not  vary.  We  have  designed  a  low-complexity  approach  for  model  order 
variation  detection  based  on  the  trace  of  the  covariance  matrix  of  the  received  signal.  Simulation 
results  demonstrate  that  the  model  order  variation  detection  approach  outperforms  those  based  on 
sliding  widow  SVD  or  MDL. 

As  multidimensional  frequency  estimation  plays  an  important  role  in  collision  resolution  when 
multiple  FH  networks  coexist,  we  have  also  proposed  an  eigenvector-based  algorithm  for  N-D 
frequency  estimation  from  multiple  data  snapshots.  We  have  analytically  quantified  the  identifi- 
ability  (in  the  noiseless  case)  and  the  performance  (in  the  noisy  case)  of  the  proposed  algorithm. 
It  is  shown  that  our  algorithm  offers  the  most  relaxed  statistical  identifiability  bound  to  date.  It 
remains  operational  when  there  exist  identical  frequencies  in  one  or  more  dimensions,  due  to  the 
adoption  of  weighting  factors.  Furthermore,  a  low-complexity  (one  iteration)  approach  is  devel¬ 
oped  to  optimize  the  weighting  factors  by  minimizing  MSEs  of  frequency  estimates.  Simulation 
results  show  that  the  proposed  algorithm  offers  better  or  competitive  performance  when  compared 
to  existing  algebraic  approaches  for  Ar-D  frequency  estimation,  but  at  lower  complexity  since  fre¬ 
quency  estimates  are  automatically  paired  without  multiple  iterations  of  joint  diagonalization  (as 
in  the  Unitary  ESPRIT)  or  requiring  an  extra  pairing  step  (as  in  MEMP  or  RARE).  It  is  shown  that 
the  optimized  weighting  factors  significantly  outperform  randomly  chosen  weighting  factors. 
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Figure  9:  (a)  Comparison  of  different  algorithms  for  3-D  frequency  estimation  from  single  snap¬ 
shot;  (b)  Comparison  of  optimized  a.  and  randomly  chosen  ra. 
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